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DIGITAL SPECTRAL ANALYSIS 


by 

Anthony J. Villasenor 
Goddard Space Flight Center 


INTRODUCTION 

Among the functions of the Test and Evaluation program at Goddard Space Flight Center are 
the following: 

1. To evaluate the structural properties of spacecraft and spacecraft components by using 
spectral statistical analysis techniques to analyze vibration test data, 

2. To simulate the spacecraft's solar environment as accurately as possible by using radia- 
tion sources with spectra that are matched with the true solar spectrum as determined by 
interferometric and spectroscopic techniques. 

Test data from these two areas, structures and solar simulation, are sent to a digital computer 
facility for spectral analysis. It was natural, therefore, that the requirements from these areas 
led to this investigation of methods for efficient spectral analysis on a digital computer. 

This paper presents some mathematical considerations of spectral analysis and some FORTRAN 
computer programs that use the Fast Fourier Transform algorithm of Cooley and Tukey. These 
programs compute Fourier amplitude and phase spectra, cross-power spectra, auto- and cross- 
correlation, and filtered spectra, as well as some other frequency domain functions. 


SAMPLING CONSIDERATIONS 

Let f(t) be a continuous function of t . The first step in analyzing f(t) by digital computer 
is to sample it and obtain a finite set of discrete points X. ( i = 1, 2, 3 . . . , n ), which will repre- 
sent the function f (t) in the computer. The process of sampling revolves around the Sampling 
Theorem (Reference 1), which states that the rate of sampling should be at least twice the maxi- 
mum frequency contained in the data. Twice the maximum frequency is actually a bare minimum, 
and experience has shown that a sampling rate of five times the maximum frequency is usually 
adequate for most applications, where the value of the maximum frequency is not precisely known. 
Insufficiently sampled signals produce spectra containing "aliased" frequencies. For example, 
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Figure 1 depicts the Fourier amplitude spec- 
trum of the curve 

f(t) = sin (277450 t /51 2). 

The spike occurs at f = 62 Hertz in the 
spectrum. The Aliasing occurs because 
the sampling rate of 512 Hertz was not enough 
to detect the true 450-Hertz signal and in- 
stead treated it as a 62-Hertz signal. We 
should have used a rate at least 2 x 450 or 
900 samples per second. The rate of 512 
samples per second gives rise to a maximum 
resolvable (cutoff) frequency f c of 256 Hertz. 
The particular value of 62 comes from f = 

2f c - 450 = 512 - 450 = 62. This relation is due to the fact that the amplitude of any frequency f in 
a signal sampled at 2f c samples per second equals the amplitude of the aliased frequency 2Nf c ± f , 
where N is the number of points in the digitized sample. Thus for a time t = f /2 


62 


256 

FREQUENCY (Hertz) 


I 

450 512 


Figure 1 -Sinusoid of 450 Hertz sampled 512 times 
per second. 


COS(27T ft) = cos 
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c 

27 
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In other words, the Fourier coefficient of the frequency f is the same as the Fourier coefficient of 
the frequency 2N f c ± f for the sampling rate of 2 f c samples per second. 


DIGITAL SPECTRAL FUNCTIONS 

Having sampled the function f(t) , we can perform various general mathematical operations 
on the data. 

The Fourier Amplitude Spectrum describes the frequency content of the data in the form of 
complex amplitude versus frequency. The term "frequency" may be misleading because it refers 
not to frequencies of our original function f(t) but rather to the number of times a particular 
sinusoid occurs within our finite sampled data. In fact, when working with a computer we are 
dealing with a Fourier Series expansion around the data even though we intend it to be a Fourier 
Integral expansion of f(t) . 

The amplitude is given by 

i \ + i I J • 


where 




-27Tik( j -1) 
e 
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The Fourier Phase Spectrum describes the phase content of our data in the form of phase 
angle (0 to 360 degrees or -180 to +180 degrees) versus frequency. The phase normally indicates 
changes of state of the f (t) generator. In the case of interferometer data, an abrupt change in 
phase means that the light source has changed its temperature relative to the detector. In vibra- 
tion data, an abrupt phase change indicates a resonance. The phase is computed as 


9 k = arc tan (I k /\) . 

The Power Spectrum is a description of the relative power of f(t) as a function of frequency. 
The Power Spectrum P(w) is the square of the Fourier Amplitude Spectrum, 


P k = R k 2 +I k 2 = (*k -^k) 

(Reference 2). The adjective "power" comes from an electrical analogy. The total energy E of a 
circuit is given by 



00 


where P is the power of the circuit in watts. But P = i 2 R , so that 


E = 



i 2 Rd t. 


If we set R to be a unit resistance, and the current i to be a time- varying function x(t) , then 
by Parseval's equality, 


E ^ 


f 

* -a 


|x(t) I 2 dt 



| F(w) | 2 dw, 


where F(w) is the Fourier amplitude spectrum of the current x(t) . The quantity |F(w)| 2 is thus 
called the power spectral density or the power spectrum. 

When x(t) is a finite, discrete function sampled at time intervals dt, we have a slightly dif- 
ferent mathematical formulation. In this case, we can speak of the average power P a over a finite 
time interval T of the current x(t) . The formulation is 
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x(t) 2 dt - 


E/T = (1/T) 


r 


If N samples are taken of the signal x( t), then T = Ndt and 


N 

p a = (l/N) 2^ x? ■ 
j = l 


We can now cite ParsevaPs equality, which holds that 



N |AJ*. 

kr 1 


Substituting this in the above, 


P a= J2 |Ak|2 ' 


Hence the average power spectral density is \\\ 2 ■ 

The power spectrum is used, for example, in acoustics to determine the system gain factor 
by finding the ratio of the output-power to the input-power spectrum. The energy or total power 
is also of frequent interest. 

The Auto-correlation function describes the degree of interaction that a function has with itself 
at various intervals along the time or distance axis. Auto-correlation for finite, discrete data is 
presented in the form of correlation coefficient versus lag number. The correlation coefficient 
at any lag varies from 1, if the function is exactly duplicated at this lag, to zero, if the function is 
completely uncorrelated with itself at this lag, and to -1 if the function cancels itself out at this 
lag. For example, any two points of a straight line are always exactly identical, hence the auto- 
correlation is everywhere 1. A sinusoid is in phase with itself once every cycle, and 180 degrees 
out of phase once every cycle, so the auto- correlation is a cosine curve with amplitude 1. Figure 
2 shows examples of auto- correlations and power spectra. 
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For a finite discrete data set x. (i = 1, 
2, 3, . . . , N ), sampled at equispaced inter- 
vals, the auto- correlation is computed as a 
function of "lag” number. For lag = 0, the 
data "comb" of N points is multiplied by it- 
self;for lag = 1, the data comb is shifted over 
one point and N - 1 multiplications are per- 
formed; for lag = 2, the comb is shifted 
again for N - 2 multiplications. The process 
continues for as many lags as desired, usually 
up to 10 percent of N. 

The auto- correlation function for finite 
discrete data is given by 


N-L 

R(L) ” R(0) N^L Xj Xj+L ’ 

j = i 


where L = lag number. Note the normalizing 
factor 1/R(0) without which this would be called an auto- covariance function. The background in 
statistics - variance, correlation coefficient - is plain. In fact, R(0) is the sample estimate of 
the true mean-square value in the data, and is related to the sample variance s 2 by R (0) = 

[(n - 1)/n]]s 2 (Reference 1). 

An important relationship exists between the auto-correlation function R (L) and the power 
spectral density P (w): they are integral transform pairs. That is, 
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Figure 2-Examples of autocorrelation and power 
spectrum functions. 


and 



P(w) e 27riw dw, 



R (L) e' 27riL dL. 


Proof of this relationship is found in the literature under the heading "Wiener- Khintchine 
Theorem" (Reference 2). This relationship is important because the usual procedure for digitally 
computing power spectra is to compute first the autocovariance and then transform it to the fre- 
quency domain to obtain a first estimate of the power spectrum. The result is only an approximation 
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because the auto- covariance was computed for a finite number of lags, whereas the true 
auto- covariance is an infinite function. 

The Cross-correlation indicates the degree of relationship between two data samples in the 
form of correlation coefficient versus lag number. Again, complete correlation is 1, complete 
uncorrelation is 0, and anti- correlation is -1. Cross -correlation differs from auto- correlation, 
being defined as an imaginary quantity. The first step in this definition is to form R^ (L) and 
R y x( L )>* where 


L = 1 ag number, 


■^xy r - — / KT T / + L , 

Vr x (0) Vr (0) 


1 N_L 

nTI ZZ Xn y " +I 


R (h) = 

yx YrJo) YT( 0) 


N-L 


N-L ZZ y " X " + L ' 


Then the even and odd parts of the cross-correlation function are given by 


E xy (L) = | [R xy (L) + R yx (L)], 

O xy (L) = \ tR xy (L) " R y x <L)] - 

The Cross Spectral Density describes the frequency dependence between two signals. If two 
signals are mutually independent and uncorrelated, the cross- spectrum is zero, but if they are 
mutually dependent, the cross- spectrum is not zero and may indicate a frequency resonance for 
the device connecting the two signals. The cross-spectrum of two signals x(t) and y(t) is a 
complex number with a real part (the co- spectrum) and an imaginary part (the quadrature spectrum) 
The co-spectrum c depicts the in-phase relationship of two signals, while the quadrature spec- 
trum Q xy depicts the out-of-phase relationship. The co-spectrum is traditionally computed as a 
cosine transform of the even part of a cross-correlation, and the quadrature spectrum as a sine 
transform of the odd part. 


CO-SPECTRUM: 


C (k) 

xy v y 


— ■i 


(0) + 2 


L. 

L 


irt k 

(-&) cos — 

L m 


+ (-l) k E(L M ) 
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QUADRATURE SPECTRUM: Q xy (k) = * J ^ 0* y (-C) sin 

N L fe 


TT^k 

1 «" 


CROSS POWER SPECTRUM: 


S = C - i Q , 

Xy xy X xy 


where the cross-correlation is + i Q^. However, it is also possible to first compute the 
cross- spectral density of two signals by direct Fourier analysis, obtain F x (w) and F y (w), then 
combine these two series to get the cross- spectral density 


S xy O) = F x (w) F y (w), 

where F y ( w ) denotes the complex conjugate of F y (w) (Reference 3). Then, by applying the Wiener- 
Khintchine theorem, we immediately get the cross-correlation function as the inverse trans- 
form of . 


FAST FOURIER TRANSFORM ALGORITHM 

Since these spectral functions involve Fourier analysis computations, any algorithm that 
offers an efficient means of computing spectra on a digital computer is much to be desired. There 
are several sources of such algorithms in the literature. Two of these appear in a book (Reference 
3). One algorithm by Goertzel makes a complex Fourier analysis of real data, using a simple 
matrix manipulation for obtaining sines and cosines. The other algorithm, by Southworth, first 
computes the auto- correlation for some percentage of the data, then performs a cosine Fourier 
transform of the result in order to obtain the power spectrum. This is the method proposed by 
Blackman and Tukey (Reference 4). 

Both of these methods have drawbacks. GoertzeTs program is exceptionally efficient, but 
still consumes too much computer time. Southworth* s program loses phase information and 
can compute only a fraction of the total frequency content of a given data set. 

More recently, James Cooley and John Tukey (Reference 5) introduced an algorithm that is 
hundreds of times faster than either of these methods. Called the Fast Fourier Transform, this 
algorithm analyzes or synthesizes a complex array of points whose size is (or can be built up 
with zeros to give) a power of 2. That is, their algorithm is set up to analyze 2 m points, where 
m is any positive integer. If the data set contains 1000 points, for example, it should be extended 
to 1024 = 2 10 points by adding 24 zero- valued points. The Cooley-Tukey algorithm is based on the 
relationship existing between two analyses of N points each and one analysis of 2N points. That is, 
if we desire a Fourier analysis of the 2N points x. ( j = 1, 2, 3, . . . , 2N ), it is faster to make one 
analysis of the odd-index points (of which there are N ) and another of the even- indexed points, and 
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then merge the two results to obtain a single analysis covering all 2N points. Since this relation 
holds between one 2N series and two N- point series, then certainly each of the N-points analyses 
can in turn be obtained by pairs of N/2-point analyses, and so on. This successive halving accounts 
for the base 2. Cooley and Tukey have estimated that this method if N/log 2 N times faster than con- 
ventional methods. Note that the errors associated with digital computation are also reduced by 
this factor of N/log 2 N. 

Thus, the plan in the Fast Fourier Transform is to generate a 2N -Point analysis from two N- 
point analyses. This procedure is straightforward. Given 2N complex data points Xj (j = 1, 2, 

3 . . . , 2N), we can split these points into an "odd" array indexed 1, 3, 5 . . . , 2N - 1, and an 
"even" array indexed 2, 4, 6 , 2N. Let sequences of A k 's and AA k 's represent the spectra 

of these arrays, respectively. Then, if we consider separately the odd and even indexed points 
in the expression for the "total" representation, by substituting 2j - 1 and 2j in place of j , we 
will eventually end with a relationship between the A k *s, the AA k r s and all the C k T s which are the 
Fourier coefficients of the 2N points. The representations would be as follows: 

Odd indexed points: 


X 


2 j _ 1 


k=l 


e 27ri(j-l)(k-l)/N ? 


Even indexed points: 


X 


2 j 


N 




e 27T i( j- l)(k-l)/N t 


where j = 1, 2, 3, . . . , N . In a 2N-point analysis, the series would be 


X. 

j 


E ^ 

k = l 


e 2 7ri( j -l)(k-l)/2N_ 


(j=l, 2,3 


, 2N) 


(The j and k indices range from 1 to N rather than the usual 0 to N - 1 to simplify programming 
notation.) If now we consider the even and odd indexed points in this last summation, 
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Odd: 


2 j 


k= 1 


e 2 77 i(2j-2)(k-l)/2N 


2N 



e 2TTi (j-l)(k - 1 )/N 


Even: 


2j 



k = i 


e 277 i (2j-l)(k-l)/2N 5 


■E 

k = 1 


c 2 tU( j=l)(k-l)/N 

H e 


e 277 i(k - 1) '2N 


Replace the index k by k' = k - N in the second half of each sum above (i.e. for k = N + l to k = 2N). 


Odd: 


2 j 


IN 

■>* E 


q 2tt i( j-l)(k-l)/N 


IN 

L Ck ' +N 


a 2?ri( j-l)(k'+N-l)/N 



e 27Ti( j_l)(k-l) N + 


E 


k ' = 1 


c e 27Ti( j-l)(k'-l)/N 

k ' + N 


since e 27T i( j *!) N/N s i. 


Even: 


■E* 

k = 1 


e 2wi(j-l)(k-l)/N 


e 2rr i(k-l)/2N 


nV 


+ 



k * - 1 


e 27ri(j-l)(k-l)/N 


e 2n i(N+k' - 1)/2N 
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k = 1 


e 27T(j-l)(k-l)/N 


e 2?r i(k-l)/2N 


N 

r 27ri(j -l)(k'-l)/N fjl 77 _ i(k # -l)/2N 

S' + N e 

k'= 1 



since 


e 27Ti (N/2N) = e" 1 = - 1. 


In these sums, k and k' are independent indices, so we can collect terms to get 


Odd: 


N 

V-.-Z 

k=l 


^k + ^k +N ^ e 


2tt i( j - 1 ) ( k - 1 )/N 


Even: 


X 2i 


= £ (C k -C k+N )e 2 - i( - 1 ^ k - 1 ^ N 


^2 77 i ( k - 1 )/2N 


Comparing this result with the definition of A. and AA we have 


A k - + C k+N ’ 


AA k ~ (^Tc " C k+N^ e 


2tt i (k-l)/2N 


Solving this system of equations finally gives 


q,* = \ 


where k = 1, 2, 3 . . . , N. 
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These equations show that it takes N complex multiplications to go from the two N-point analyses 
to a 2 N-point analysis. To find the total time needed to make an N-point analysis by this method, 
let M n equal the number of multiplications needed for each of our N-point analyses. Let M 2N equal 
the total number of multiplications needed to end with the 2 N-point analysis. Then obviously 

M 2N 3 2M n + N. 

Since N = 2 m , m successive doublings would give us the required N-point analysis - that is to say, 
for the above iterative equation the initial condition is N 0 = 1. Solving the equation 


M 2Ni = 2[ y + N i 

(i = 0, 1, 2 . . . , m- 1) we get M n = (n/ 2) log 2 N . To compensate for various computer overhead 
operations, this can be overestimated as Nlog 2 N operations, as opposed to N 2 operations by 
Goertzel’s method. Hence the saving in computer time using Cooley’s method is N/log 2 N times 
the time by the old methods. 

For a Control Data Corporation (CDC) 3100 computer with software floating point, Goertzel’s 
algorithm took 45 minutes to Fourier-analyze 2048 points. It took Cooley’s method 18.3 seconds — 
an enormous saving. It makes possible real-time digital spectral analysis. See Table 1 for 
comparisons of execution times. 


Table 1 


Table of Comparison Times in Seconds 
(CDC 3100 with Software Floating Point) 




N 

size of array 




256 

7 

512 

1024 

2048 

4096 

8192 

Goertzel 

42.9 

169.8 

675.7 

2699.5 



Harmon 



14.0 




Fourier 

1.7 

3.8 

8.4 

18.3 

31.9 


Aliasing 

2.7 

5.5 

11.5 

24.2 



Double 



18.9 


62.8 

114.4 


COMPUTER PROGRAMS 
Subroutine FOURIER 

The original subroutine HARMON written by Dr. Cooley to carry out these operations is some- 
what limited because it requires a complex data array as input. Most numerical data are real, not 
complex. To solve the majority of problems, the subroutine FOURIER was written to take as 
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input a real N -point array and to return (n /2) + 1 complex amplitudes. Subroutine FOURIER 
accomplishes this by the following procedure: First, it stores the odd-indexed points of the given 
N -point array in the real parts of n/ 2 complex numbers, and the even-indexed points in the 
imaginary parts. <Next, it makes a call to subroutine HARMON with n/ 2 complex points as input. 
Finally, it extracts the odd and even pairs of amplitudes from the result and merges these to form 
the (n/ 2) + 1 complex amplitudes of the original N points. The number of complex amplitudes is 
(N/2) + 1 because the Fourier coefficients for real data satisfy: 

\ =^- k + 2 ( k = 1* 2, 3 . . . N, N + 1) . 

That is, the real and imaginary parts are respectively symmetric and antisymmetric about the 
element k = (n/ 2) + 1. This relation yields the following identities: 

Real part of = average of input data points 
Imaginary part of A 1 = 0 
Imaginary part of A (N/2) + t = 0 
\ = A n ^ x (periodicity). 

Thus, for real analysis we need only be concerned with the (n/ 2) + 1 complex amplitudes A v 
A 2 , . . . , A n/2 A (N/2)+1 . At this point, the procedure for accomplishing the real analysis is 
straightforward. Given a complex array x k = x k + i x k , where the superscript represent the odd 
and even indexed real data points, respectively, we call HARMON to get the complex amplitudes 
C k (k = 1, 2, 3 . . . , n). Thus any original complex point x k has the series representation 


N 



( 1 ) 


At this point we can also introduce the Fourier coefficients A and AA. ( j = 1, 2, 3 . . . , n) repre- 
senting respectively, the odd- and even-indexed series coefficients. We can write 


x° 


E 

i =1 


A e 27ri(j-l)(k-l)/N 
3 


( 2 ) 


and 


x 


e 

k 


E 

i=l 


AA. 

j 


e 277 i(j-l)(k-l)/N 


(3) 
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Substituting Equations 2 and 3 in Equation 1 gives 


X° + iXt - ^ 1 1 in \ (j-l)(k-l)/N 


± iX k = / (A. ± i AA.) e 


j = i 


Taking the complex conjugate in Equation 1 gives 


(4) 


i =1 


c. 


e - 277 i(j-l)(k-l)/N 



2 77 i (" £ - l)(k-l)/N 
e 


(5) 


where in the second summation the index j has been replaced by l = N - j + 2. That is, starting 
with 


N 

£ e ~2rr i *(j-l)(k-l)/N 
J 

J =1 


and letting j = N - l + 2, immediately gives the following identities: 

C j = 2 

<2w i(j-l)(k-l)/N — 2 , rr i ( N — 't + 1 ) ( k — 1 ) /N 2ffi(^ - l)(k-l)/N 

e = e - e 


j = 1 becomes = N + 1 
j = N becomes - 2 




j= i 


becomes 





( N+l ) + 2 


= ^N+l ~ C 1 


so that (C) is demonstrated. 



Combining Equation 1 and 4 gives 


C. = A. 

j j 


+ iAA. 


Combining Equations 4 and 5 gives 




- j + 2 


A. 


iAA. 

j 


Thus, 


A i = 7( c i 

AA i =-~ (C i + 


Now we have the odd and even amplitudes, which can be recombined to give 


C. 

j 


— (A. + AA, e 2 
2 1 1 


i C j -t)/N 


)• 


C. - — (A. - AA . e ^i(i-D/N 1 A aa. e' 2,Ti(i - 1)/N ), 
j 2 J 1 2 V 1 } J 


which is the desired result. 

Subroutine ALIASING 

Another routine written around the original Cooley program is subroutine ALIASING. Its 
purpose is to determine, from N given real points, whether the sampling rate is sufficient. It also 
computes ParsevaFs formula to determine the magnitude of calculation errors attributable to 
digital roundoff. The routine returns with (n/ 2) + 1 complex amplitudes, just as if the routine 
FOURIER had been called instead. 

The method used in ALIASING is similar to the one in FOURIER: An N-point array is treated 
as an (N/2)-point complex array, which is Fourier-analyzed with HARMON. The resulting spectrum 
is split up into two spectra A k and representing, respectively, the frequency components of the 
odd and even indexed points of the given array. Each of these spectra represents a sampling 
density of N/2 points; that is, half the actual density. Then both the odd and the even spectra 
are merged to yield the desired final spectrum C k . The A k f s and AA k t s represent the same sampling 
density at alternating points, so that their difference (between odd- indexed and even-indexed points) 
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should be zero for each k . If the sampling rate is sufficient, the differences between A k and c k 
and between AA k and C k will also be close to zero. In addition, all the octaves c N _ k should be zero, 
i/ any or all of these quantities are larger than the difference between A k and AA k , then the sampling 
density is insufficient or, at best, questionable. 

The quantities returned by subroutine ALIASING are computed by the following formulas: 


Odd- even: 


N/4 



l\l - l A \l 


Odd-both: 


N/4 

L 


k = l 


l\l 


1^1 


2 


Even-both: 


N/4 

L I |a \i - 1^ 1 

k=i 


Octave: 


N/4 



S ^/2 


-k + 2 


2 


As an example, suppose that we generate a unit sinusoid of 450 Hertz, and sample it at 1024 
points per second for 1 second. We know that all C k r s will be zero except c 450 , which will equal 
1. A call to ALIASING will yield zero amplitudes for all A T s and AA's except for k = 62, when both 
will equal 1. The frequency of 62 Hertz is the aliased frequency corresponding to 450 Hertz for 
a sampling density of 512 points per-second. However, c 62 will be zero, while c 450 will be 1. This 
procedure would show the difference between the A r s and AA T s as zero, but obviously not zero for 
the other differences. In order to capture the 450- Hertz sinusoid, we would have to increase the 
sampling rate, say by a factor of 2. 

Parseval's theorem simply says that the norm of the input data array should be equal to the 
norm of the output amplitude array. This is a property of linear transformations used in computer 
applications to determine the extent of error resulting from digital truncation. The statement of 
ParsevaTs equality is 
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N 

E 


X. 


N 

E 

k= 1 


\ 


ALIASING deals with real data, however, so that 

\ ~ ^N-k+2‘ 

Also, the amplitude of all components but \ and A (N/2) + 1 are doubled because the routine computes 
a symmetric two-sided spectrum, but only positive frequencies are of interest. Hence Parseval's 
equality is computed as 


D* 

i = l 



|Aj 2 + 2 (A 2 


+ A (N/2) + 


Subroutine DOUBLE 

It sometimes happens that the sampling rate yields too many points for a particular computer. 
A Control Data Corporation (CDC) 3100 with 16,000 words of memory and standard FORTRAN can 
handle at most 4096 points, using the Fast Fourier Transform. The routine called DOUBLE was 
written to double any computers capacity for Fourier analysis; it requires three scratch tapes. 
The user puts his oversized array on tape, makes a call to DOUBLE, then reads in that same 
tape to get his complex amplitudes. The user of course does not read and write his oversized 
N -point array in a single command; instead, he must divide his array into four records of n/ 4 
points each. DOUBLE does the same processing as FOURIER except for input/output (I/O) calls 
and increased capacity. It should be mentioned that the complex amplitude A^ N/2) + 1 is not returned 
from DOUBLE. Also, if only frequency amplitudes are required and phase information is dis- 
carded, then the user can save I/O usage and computer time by modifying DOUBLE to return just 
the n/2 amplitudes. 

Subroutine POWER 

Another subroutine power takes as input data two real arrays X and Y and returns in x with 
the cross-power spectrum 


p xy (w) = 3 (X) 3 (Y) . 
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The power spectral density of a single array X may be computed by power as 


P xx (w) =3(X) 3(X) = |3(X)| 2 , 

but this would not be as efficient as computing the squared modulus of the complex spectrum re- 
turned by FOURIER. 

Subroutine COVAR 

The routine COVAR uses CROSS to yield the cross- covariance of two given arrays. The co- 
variance is computed as the inverse Fourier transform of the cross-power spectrum. COVAR 
first computes the cross-power spectrum by calling CROSS, then performs an inverse Fourier 
transform on this array to get the cross-covariance. For input arrays x andY, the cross covariance 
is returned in Y and the cross power spectrum in X . To obtain the correlation function, normalize 
the covariance values by the "zero lag" value of the covariance - i.e., the value Y x . 

Subroutines HANNING and HAMMING 

The subroutine HANNING takes as input a real array A of frequency amplitudes. The Hanning 
smoothing function is applied to the input, and the result is stored in array A, thereby replacing 
the original contents of A. The Hanning function is defined by 

A' (1) = A(l) + A(2), 

A'(N) = A(N - 1) + A(N), 

A' (i ) = . 5 A(i - 1) + A(i) + . 5 A(i +1) 

(i r 2, 3, 4. . . N _ 1). 

The coefficients are customarily given as (0.25, 0.5, 0.25), but we are dealing here with one-sided 
spectra so that a factor of 2 is involved. 

The subroutine HAMMING applys the Hamming smoothing function to a given real array A . It 
is defined by 


a; 3 1.08 Aj + 0.92 A 2 , 

K = 1-08 A* + 0.92 Vr 
A'. = 0.46 A._ x + 1.08 A. + 0.46 A. + x 
(i = 2, 3, ... N - 1). 

Note that the factor 2 appears again in the coefficients, which are usually given as (0.23, 0.54, 0.23). 
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The usefulness of these little subroutines for analyzing data with discrete frequencies should 
not be underrated; they compensate for the finite discrete nature of Fourier analysis when per- 
formed on a digital computer. These simple digital filters perform two functions: they reduce 
the resolution of spectral lines, and they increase the accuracy of relative height of spectral lines. 
For continuous, smooth spectra, this is not too important, but spectra with sharp spectral lines can 
be grossly misinterpreted if these two functions are not performed. In spite of the reduction in 
resolution, some filter should be used in digital Fourier analysis. No filter actually corresponds 
to using asinx/x filter, which has undesirable high sidelobes that can lead to gross errors of com- 
puted results. 


AVERAGING FOR BETTER RESOLUTION 

Often, a sample size of n points does not adequately describe the frequency content of the data, 
probably twice as many points being needed. At this point it is often suggested that, rather than 
resample the data at a higher sampling rate, we should just perform linear interpolations between 
the given points to get the required point density - say, 2N points. The reply is that a linear 
interpolation - or rather, a linear transformation - adds no new information to the data, hence no 
new information to the spectrum. The proof is as follows: 

Consider the N-point series representation 


IN 

■.-z> 


e 2wi ( j -l)(k-l)/N 


Suppose that we define a new sequence of points XX. located between the x. , such that 

XX. = a. X. + b. X. , 
j j j j j+i 

and such that the corresponding Fourier coefficients AA . are represented by 


XX. = £ A \ 


In order to get the AA^s in terms of the A k 's, we set 
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X 


j+1 


L 

k=l 


\ 


e 2 7T i j(k-l)/N 


■£ 

k = l 


\ 


e 27Ti( ] - l)(k-l)/N 


g 2tt i(k-X)/N 


Hence, 


XX. = a. 
j j 


IN 


2ni(i -l)(k-l)/N + b 


IN 

2 > 


2rr i( j-l)(k-l)/N 27ri(k-l)/N 


k=l 


■£ 

k - 1 


( a i \ + b j \ 


q 2tt i(k-l)/N^ 


e 2ffi(j-l)(k-l)/N 


or 


AA, = a, \ + b. A, e 2 ^ k ' 1)/N 

This equation can easily be made independent of j by assuming a uniform linear interpolation such 
that a = a. and b = b. # Then we have 

A\ = aA, + b/V « ? " l(k - 1>/N 

* , , 2 7Ti(k-l)/N x 

= \ (a + b e )■ 

Now if we merge these points X. and XX. , we form a 2N-point series X' , such that our original 
points are the odd indexed x' and the linearly interpolated points are the even indexed X'. . Using 
the equalities obtained earlier for merging series, we get, for k = 1, 2, 3, . . . , N + 1. 


q =1 (\ + aa. 


e -2iri(k-l)/2N^ 


2 


[\ + \ (a + b e 277 


(k-l)/Nj e 


- 2n i(k-l)/2N 


] 


1 

2 


\ [ 1 + (a + b e 


2 n i<k- 1)/N 


-2-tt i(k— 1)/2N ■ 
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Also, 


C N + k = \ \ [l-fa+be 2 '^") e -2w i(k-l)/2N] 

This shows that the new series, with alternating points being linearly interpolated values, is merely 
a modification of the old series by a complex exponential whose period is 2N points. The old fre- 
quency components will be modified and no new ones will appear. This effect is easily seen if we 
set a = b = 1/2, for then 




1 +_ (1 + e 27Ti(k-l)/N) e -27Ti(k-l)/2N 
2 




e - 2 77 i(k-l)/2N + e 2 7T i(k-X)/ 2N 
2 




+ cos 


277(k - 1) 
2N 


and 




+k 


1 _ 

2 




cos 


277(k - 1) 
2N 


The curve f(t) defined by 




> 


(k= 1, 2, 3. . . , 
N+l) 


J 



(b) SPECTRUM (typical) 


r 


1 

2 


1 + cos 


2 "77 (k - 1) 

2N 


for first half 
of spectrum 


f(t) < 


2 


cos 


277 (k-1) 

2N 


for second hal f 
of spectrum 



(c) SPECTRUM WITH CURVE f (t) 

Figure 3— Effect of averaging neighboring points 
in a data sample. 


(k = 1, 2, 3... , N+ 1) 

has the form shown in Figure 3(a). A typical 
spectrum of N points is shown in Figure 3(b). 
The result of linearly interpolating between 
points to get 2N points is shown in Figure 3 (c) . 
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DIGITAL FILTERS 


When an infinite and continuous function h(t) is sampled for digital processing, the actual function 
considered by the computer is d(t) = g(t)h(t) , where g(t) is defined by 

M 

g(t) = 8 (t - j A t). 

j=-M 

It equals 1 when t = j At and is otherwise zero. This finite sequence of N, N = 2M + 1 , equispaced 
unit spikes, corresponding to N sampled points, is called a finite Dirac comb. Function g(t) is a 
specimen of the so-called "lag windows" that modify h(t) at different intervals or time lags. 

When a Fourier transform is applied to d(t), the result is the same as convolving the Fourier 
transforms of g(t) and h(t). That is, if G(w) and H(w) are the Fourier transforms of g(t) and h(t), 
respectively, then the Fourier transform of the product g(t)h(t)is the convolution G*H, defined as 

G*HO) - f G(f)H(w- f)clf - H*G(w). 


To see this, we compute 



g(t)h(t)e- 2ffi wt dt 


r° 

-CO 


h(t) 

G( f ) e 2 ’ f ‘ 

1 df 

J~co 

j -00 

_ 


e -2-i wt dt 



df 


Thus, instead of multiplying two functions before transformation, we could equivalently convolve 
their spectra represented by 


G(w) 



g(t) e- 


wt dt 
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and 


H(w) = I h(t) e- i5,t dt. 

•'-00 

The convolution G* H(w) represents an estimate of the true frequency content H(f) as seen through 
a spectral window of variable transmission g( w - f). This spectral window is referred to as a 
digital filter when applied to digital data; in that case, the discrete convolution takes the form 

M 

G*H(n) = J2 G(j) H(n - j), 
i 

which has the effect of spreading the "true" spectral amplitude H(n) by M weights on either side of 
n. The H(n) are called "raw spectral estimates"; the G*H(w) are "refined (or smoothed) spectral 
estimates". In most cases G(w) has the form of a (sin L)/L function, which has a maximum peak at 
L = 0 (equal to 1} and an infinite number of relative maxima oscillating in sign. The peak at L = 0 
is called the main lobe, and the relative maxima are called side lobes. 

The choice of a particular digital filter is based on the following considerations: 

1. g(t) - 0 for all t > | T | 

^ 0 otherwise 


2. The main lobe of G(w) should be concentrated near L = 0 to reduce the resolution or band- 
width represented by a particular frequency at w. This means that g(t) should be flat with sharp 
cutoffs at ±T, such as a square wave. 

3. The side lobes should be as small as possible to reduce the contribution of other frequen- 
cies at w - f to the final amplitude at w. This means that g(t) should be smooth and gradually 
trailing off the 0 at t = ± T. 

The only way to satisfy all these requirements is to effect some sort of compromise deter- 
mined by experiment. 

The HAMMING and HANNING subroutines described earlier are digital filters that smooth 
the spectrum and give discrete amplitudes a better relative representation at the cost of individual 
peak resolution. It would be worthwhile to see directly the relationship between lag windows and 
spectral windows for the popular Hanning function. The Hanning lag window is defined as 
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1 


/ n 1 I i 77 t 

g(t) = *2 l 1 + cos T" 


= 0 


t < T 


t > T 


Figure 4 is the graph of the function. 

To obtain the Fourier spectrum of g(t) , we 
compute 


G(f) 


fCO 
" -00 


g(t) e ~ 2ni ft dt 


if 

2T J. 



dt + 


if g(t)e- 2 ’ 7ift dt 
J_T 

1 f T TTt _ 2 77 i 

I cos e 

2 T J _ T 


TIME, t (sec) 

Figure 4— Hanning lag window function in 
time domain. 


dt 


/ T 

(CO 


s277ft-i sin 277ft)dt + 


-l f 

2T i. 


77 t 

cos (cos 277ft- i sin 27rft) dt 

T 


sin 2 77 f T 1 
27 Tf + ~2 


in 2 77 ({ - 


T i sin 277 ( f + ^ ]T 


J. 

2 T 


277 f - 

2 T 


2 77 f + 


2 T 


Thus the spectral form of the Hanning lag window is the sum of a central term at frequency f , and 
two other (sin x )/x terms displaced on either side of f by 1/2 T. The difference 


f t J_\ - U - JL\ - — 

2 T / \ 2 T j T 


is the digital bandwidth of the filter. For a sample of N points, the bandwidth is 2/(NAt), since 
NAt = 2 T, Figure 5 is the graph of G(f) . Substituting G(f) in the convolution formula, we have 


G*H( f ) 



H(f - l) dl 



H(f--t) 




+ 


2 


sin 277 


277 



— w 

2 T J 


d l . 


In the discrete case, a frequency f = f represents the k th harmonic in a finite sample of length 
2 T , so the f k can be expressed as f k = k/ 2 T . Hence if f = f k and^t = f ^ we get 
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NORMALIZED NORMALIZED 

AMPLITUDE AMPLITUDE 



Figure 5— Hanning lag window function in 
frequency domain. 



Figure 6— Dirac spectral window function. 




1 sinTT^-l) 1 sin7T(^+l) 

2 77 (Z - 1) 2 77 ('f' + 1) 


Integrating by parts gives 

G*H(f k ) =1 H(f k ) +I H (f k+1 ) + i Hff,.,). 

The coefficients are the triplet (0.25, 0.5, 
0.25) usually given for the Hanning spectral 
window. If we are dealing with one-sided 
spectra (i.e., 0 < f < x) rather than two-sided 
spectra (-00 < f < cc), then these coefficients 
should be doubled. 

The spectral window corresponding to 
the finite Dirac comb is 


G(f) = NAt cos 



sin 2^/k 

277k 

sin 77k/N 


77k 

N 


Figure 6 is the graph of this function. The lobes are not as small as in the Hanning window, and 
for that reason the latter is usually preferred. 

Further information on digital filters can be found in Reference 4. 


Goddard Space Flight Center 

National Aeronautics and Space Administration 
Greenbelt, Maryland, December 12, 1967 
125-23-02-13-51 
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Appendix A 

Subroutine HARMON 
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SUBROUTINE HARMON ( A , S ,M » I FS , I FERR ) 

DIMENSION A(l) .5(1) 

HARM, ONE-DIMENSIONAL BASIC FORTRAN VERSION. J.W. COOLEY HARM 001 

MODIFIED TO RUN ON CDC 3100. 

HARM 002 
HARM 008 

DOES EITHER FOURIER SYNTHES I S , I . E . , COMPUT ES COMPLEX FOURIER SERI ESHARM 009 
GIVEN A VECTOR OF N COMPLEX FOURIER AMPL I TUDES , OR , GIVEN A VECTOR HARM 010 


OF COMPLEX DATA X DOES FOURIER ANALYSIS, COMPUTING AMPLITUDES. HARM Oil 

A IS A COMPLEX VECTOR OF LENGTH N=2**M COMPLEX NOS. OR 2*N REAL HARM 012 

NUMBERS. A IS TO BE SET BY USER. HARM 013 

M IS AN INTEGER 0 • LT .M . LE • 1 3 * SET BY USER. HARM 014 

S IS A VECTOR S ( J ) = SIN { 2*P I * J/NP ), J=l,2» NP/4~1, HARM 015 

COMPUTED BY PROGRAM. HARM 016 

I FS IS A PARAMETER TO BE SET BY USER AS FOLLOWS- HARM 017 

I F S — 0 TO SET NP = 2**M AND SET UP SINE TABLE S. HARM 018 

I FS= 1 TO SET N=NP=2**M, SET UP SIN TABLE. AND DO FOURIER HARM 019 

SYNTHESIS, REPLACING THE VECTOR A BY HARM 020 

HARM 021 

X(J)= SUM OVER K = 0*N-1 OF A ( K ) *EXP ( 2*P I * I /N ) ** ( J*K ) , HARM 022 

J = 0 , N- 1 , WHERE I =SQRT ( - 1 ) HARM 023 

THE X-S ARE STORED WITH RE XU) IN CELL 2*J+1 HARM 024 

AND IM X(J) IN CELL 2*J+2 FOR J = 0 , 1 , 2 * . . . * N-l . HARM 025 

THE A-S ARE STORED IN THE SAME MANNER. HARM 026 

HARM 027 

I FS=-1 TO SET N=NP = 2**M,SET UP SIN TAB*-E, AND DO FOURIER HARM 028 

ANALYSIS, TAKING THE INPUT VECTOR A AS X AND HARM 029 

REPLACING IT BY THE A SATISFYING THE ABoVE FOURIER SERIES. HARM 030 

I FS=+2 TO DO FOURIER SYNTHESIS ONLY, WITH A PRE-COMPUTED S. HARM 031 

I FS=— 2 TO DO FOURIER ANALYSIS ONLY, WITH A PRE-COMPUTED S. HARM 032 

IFERR IS SET BY PROGRAM TO- HARM 033 

=0 IF NO ERROR DETECTED. HARM 034 

=1 IF M IS OUT OF RANGE., OR, WHEN IFS=+2*“2, THE HARM 035 

PRE-COMPUTED S TABLE IS NOT LARGE ENOUGH. HARM 036 

=-l WHEN I FS =+1,-1, MEANS ONE IS RECOMPUTING S TABLE HARM 037 

UNNECESSARILY. HARM 038 

HARM 039 

NOTE- AS STATED ABOVE, THE MAXIMUM VALUE OF M FOR THIS PROGRAM HARM 040 

ON THE IBM 7094 IS 13. ON 360 MACHINES HAVING GREATER STORAGE HARM 041 

CAPACITY, ONE SHOULD CHANGE THIS LIMIT By REPLACING 13 IN HARM 042 

STATEMENT 3 BELOW BY LOG2 N, WHERE N IS THE MAX. NO. OF HARM 043 

COMPLEX NUMBERS ONE CAN STORE IN HIGH-SPEED CORE. HARM 044 

IF THE CAPACITY OF HARM IS TO BE INCREASED, ONE MUST HARM 045 

ALSO ADD MORE DO STATEMENTS TO THE BINARY SORT ROUTINE HARM 046 

FOLLOWING STATEMENT 24 AND CHANGE THE EQUIVALENCE STATEMENTS HARM 047 

FOR THE K-S. HARM 048 

HARM 049 


DIMENSION K ( 12 ) 

EQUIVALENCE ( K ( 1 1 ) , Ki ) , ( K (10 ) ,<2 ) , < K < 9 ) , K3 ) * { K ( 8 ) , <4 ) , ( K ( 7 ) ,K5 ) 
EQUIVALENCE ( K < 6 > ,<6 ) * < K ( 5 ) * K7 ) , < K { 4 ) *Ke ) * ( K ( 3 ) ,K9 ) , ( K ( 2 ) *K 1 0 > 
EQUIVALENCE ( K ( 1 ) , K 1 1 ) , ( K ( 1 ) , N2 ) 

IF(M)2*2»3 
3 IF(M-ll) 5*5,2 
2 I FERR=1 
1 RETURN 
5 I FERR=0 
N=2**M 

I F ( IABS(IFS) - 1 ) 200,200,10 

WE ARE DOING TRANSFORM ONLY. SEE IF PRE-COMPUTED 
S TABLE IS SUFFICIENTLY LARGE 
10 I F ( N-NP ) 20,20 , 12 


HARM 055 

HARM 057 
HARM 058 
HARM 059 
HARM 060 
HARM 061 
HARM 062 
HARM 063 
HARM 064 
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n n 


12 

IFERR=1 

HARM 065 


GO TO 200 

HARM 066 


SCRAMBLE A* BY SANDE-S METHOD 

HARM 067 

20 

K(1)=2*N 

HARM 068 


DO 22 L = 2 * M 

HARM 069 

22 

K(L)=K(L-l)/2 
DO 24 L=M » 10 

HARM 070 

24 

K ( L+l ) =2 

HARM 072 


NOTE EQUIVALENCE OF KL AND KI14-L) 

HARM 073 


BINARY SORT- 

HARM 074 


I J=2 
J 1 = 2 

HARM 075 

25 

DO 30 J2 = J1 *K2 * K1 



DO 30 J3=J2»K3*K2 

HARM 078 


DO 30 J4=J3»K4*K3 

HARM 079 


DO 30 J5=J4»K5*K4 

HARM 080 


DO 30 J6 = J5 *K6 *K5 

HARM 081 


DO 30 J7 = J6 * K7 * K6 

HARM 082 


DO 30 J8=J7*K8*K7 

HARM 083 


DO 30 J9=J8,K9,K8 

HARM 084 


DO 30 Jl 0 = J9 »K1 0 , K9 
DO 30 J I =J1Q*K1 1 *K10 

HARM 085 


I F ( IJ-JI )28,30*30 

HARM 089 

28 

T = A ( I J — 1 ) 

HARM 090 


A ( IJ — 1 > = A(JI — 1 ) 

HARM 091 


A { J I -1 ) =T 

HARM 092 


T=A( IJ) 

HARM 093 


A(IJ)=AI JI ) 

HARM 094 


A ( J I )=T 

HARM 095 

30 

I J= I J+2 
J 1= Jl+2 

IF(K1-J 1)31*25*25 

HARM 096 

31 

I F ( I FS ) 32 * 2 * 36 

HARM 097 


DOING FOURIER ANALYSIS, SO DIV. BY N AND CONJUGATE. 

HARM 098 

32 

FN = FLOAT ( N ) 



DO 34 I = 1 *N 

HARM 100 


A(2*I-1 ) = A(2*I-1)/FN 

HARM 101 

34 

A(2*I ) =— A (2*1 ) /FN 

HARM 102 


SPECIAL CASE- L=1 

HARM 103 

36 

DO 40 1=1, N, 2 

HARM 104 


T = A ( 2*1 “1 ) 

HARM 105 


A (2*1—1 ) =T + A ( 2* I + 1 ) 

HARM 106 


A ( 2* I +1 ) =T-A (2*1+1 ) 

HARM 107 


T =A ( 2* I ) 

HARM 108 


A(2*I ) = T + A ( 2* I +2 ) 

HARM 109 

40 

A ( 2*1+2 ) = T - A ( 2* I +2 ) 

HARM 110 


IF(M-l) 2*1 ,50 

HARM 111 


SET FOR L=2 

HARM 112 

50 

LEXP1 =2 

HARM 113 


LEXPl=2**(L-l ) 

HARM 114 


LEXP=8 

HARM 115 


LEXP=2**(L+1 ) 

HARM 116 


NPL= 2**MT 

HARM 117 


NPL = NP* 2**-L 
DO 130 L=2*M 

HARM 118 


SPECIAL CASE- J=0 

HARM 120 


DO 80 I =2 *N2 » LEXP 

HARM 121 


11=1 + LEXP1 

HARM 122 


12=11+ LEXP 1 

HARM 123 


13 = I 2+LEXP 1 

HARM 124 


T=A( 1-1 ) 

HARM 125 


A ( I -1 ) = T +A( 12-1 ) 

HARM 126 


A ( 12-1 ) = T-A( 12-1 ) 

HARM 127 


T =A( I ) 

HARM 128 


A ( I J = T + A( 12) 

HARM 129 


A ( I 2 ) = T-A ( 12 ) 

HARM 130 


T = -A( 13) 

HARM 131 


T I = A ( I 3-1 ) 

HARM 132 


A ( I 3-1 ) = A(I1-1) - T 

HARM 133 


A ( I 3 ) = A(I1 ) - T I 

HARM 134 


A ( 1 1 — 1 ) = A ( I 1- 1 ) +T 

HARM 135 
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I 


80 Aim = A(I1 ) + TI 

lF(L-2) 120*120*90 
90 KLAST=N2-LEXP 
JJ=NPL 

DO 110 J = 4 > LEXP 1*2 

NPJJ=NT-JJ 

UR = S { NP J J ) 

U I =S { J J ) 

I LAST = J+KLAST 
DO 100 1= J * I LAST * LEXP 
I1=I+LEXP1 
12=1 1+LEXP 1 
I3=I2+LEXP1 
T=A( 12 — 1 ) *UR-A ( 12 )*UI 
T I = A ( 12-1 )*UI+A( I 2 ) *UR 
A < 12-1 ) = A ( I -1 ) -T 
A ( 1 2 )=A( I ) - TI 

A ( I — 1 ) =A ( I -1 ) +T 
A ( I ) =A ( I ) +T I 
T = -A( I 3- 1 ) *U I -A ( I 3 ) *UR 
T I =A ( I3-1)*UR-A( I 3 ) *U I 
A ( 13-1 )=A( 1 1-1 )-T 
A ( I 3 ) =A ( 1 1 ) — T I 

A ( I 1-1 )=A< 1 1-1 )+T 
100 A ( 1 1 ) =A( ID +TI 

C END OF I LOOP 

110 J J= JJ+NPL 
C END OF J LOOP 

120 LEXP1=2*LEXP1 
LEXP = 2*LEXP 
130 NPL=NPL/2 
C END OF L LOOP 

I F ( IFS) 143*2*1 

CC DOING FOURIER ANALYSIS. REPLACE A BY CONJUGATE. 

145 DO 150 I = 1 * N 
150 A ( 2* I ) =-A ( 2* I ) 

GO TO 1 
C RETURN 

C MAKE TABLE OF S ( J ) =S I N ( 2*P I * J/NP ) * J=1 *2 * . • . .NT-1 *NT=NP/4 

200 NP=N 
MP=M 
N T = N/4 
MT = M-2 

IF(MT) 260*260*205 
205 THETA=. 7853981634 
C THETA=PI/2**(L+1 ) FOR L=1 

JSTFP = NT 

C JSTEP = 2**1 MT-L+1 ) FOR L=1 

JDIF = NT/2 

C JDIF = 2**(MT-L) FOR L=1 

S( JDIF) = SIN( THETA) 

IF (MT-2)260*220*220 
220 DO 250 L = 2 * MT 

THETA = THETA/2. 

JSTEP2 = JSTEP 
JSTEP = JDIF 
JDIF = JDIF/2 
S( JDIF)=SIN(THETA) 
jci=nt-jdif 
S( JC1)=C0S( THETA) 

JLAST = NT -JSTEP 2 
IF( JLAST- JSTEP) 250*230*230 
230 DO 240 J = JSTEP * JLAST *JSTEP 
JC=NT- J 
JD= J+ JD I F 

240 S( JD)=S{ J)*S( JC1 )+S{ JDIF)*S( JC) 

250 CONTINUE 
260 IF( IFS)20*1*20 
END 


HARM 136 
HARM 137 
HARM 138 
HARM 139 
HARM 140 
HARM 141 
HARM 142 
HARM 143 
HARM 144 
HARM 145 
HARM 146 
HARM 147 
HARM 148 
HARM 149 
HARM 150 
HARM 151 
HARM 152 
HARM 153 
HARM 154 
HARM 155 
HARM 156 
HARM 157 
HARM 158 
HARM 159 
HARM 160 
HARM 161 
HARM 162 
HARM 163 
HARM 164 
HARM 165 
HARM 166 
HARM 167 

HARM 169 
HARM 170 


HARM 173 
HARM 174 
HARM 175 
HARM 176 
HARM 177 
HARM 178 
HARM 179 
HARM 180 
HARM 181 

HARM 183 
HARM 184 
HARM 185 
HARM 186 
HARM 187 
HARM 188 
HARM 189 
HARM 190 
HARM 191 
HARM 192 
HARM 193 
HARM 194 
HARM 195 
HARM 196 
HARM 197 
HARM 198 
HARM 199 
HARM 200 
HARM 201 
HARM 202 
HARM 203 
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Appendix B 

Subroutine FOURIER 


SUBROUTINE FOURIER (A»S»M*IFS) 

THIS ROUTINE PERFORMS AN ANALYSIS OF 2**M POINTS BY FIRST DOING 
AN ANALYSIS OF 2**M/2 COMPLEX POINTS AND THEN ARRANGING THE RESULTS 

ARGUMENTS 

1. A - REAL DATA ARRAY - OF DIMENSION 2**M + 2 

2. S - SIN/COS TABLE - DIMENSION 2**<M-3) 

3. M - EXPONENT OF 2 - SIZE OF REAL ARRAY 

A. I FS 1 FOR FIRST TIME* -2 THERAFTER 

DIMENSION A ( 1 ) »S ( 1 ) 

N = 2**(M-1) 

CALL HARMON I A * S * M-l * IFS* IFERR ) 

C MERGE 2 N-POINT ANALYSIS INTO I 2N-POINT ANALYSIS 

NHALF = N/2 
NTWO = N*2 + 4 

X = XO = COS( 3* 1415926536/FLOAT ( N ) ) 

Y = YO = S IN ( 3 • 1415926536/FLOAT ( N ) ) 

DO 1000 K2 = 4 * N * 2 

K1 = K.2 - 1 

N2 = NTWO - K2 

N 1 = N 2 - 1 

BK1 = A ( K 1 ) + A ( N 1 ) 

BK2 = A ( K2 ) - A I N2 ) 

BN1 = A(K2) + A(N2) 

8N2 = A I K1 ) - A I N 1 ) 

XBN1 = X*BN1 
XBN2 = X*BN2 
YBN1 = Y*BNl 
YBN2 = Y*BN2 

A ( K1 ) a .5 *(BK1 + XBN1 - YBN2 ) 

A(K2) = .5 * ( -BK2 + XBN2 + YBNl ) 

A ( N 1 ) = .5 *<BK1 - XBN1 + YBN 2 ) 

A(N2) = .5 *<BK2 + XBN2 + YBNl) 

Q = X*XO - Y*YO 

Y = Y*XO + X*YO 
1000 X = Q 

C COMPLEX ELEMENT AIN) 

A ( 2*N+1 ) = ( A ( 1 ) - A ( 2 ) ) * • 5 
A(2*N+2) = 0,0 
C COMPLEX ELEMENT A(O) 

All) = • 5* I A ( 1 ) +A ( 2 ) ) 

A I 2 ) = 0.0 

C COMPLEX ELEMENT AIN/2) 

C AIN+l) = AlN+l) 

C AIN+2) = AIN+2) 

RETURN 

END 
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Appendix C 

Subroutine ALIASING 


subroutine ALIASING ( a*s*m* DATASUM* AMPSuM * oddevens. oddboths* 

1 BOThSUM* OCTAVES) 

this program compares the results of the given sampling with 

ARRAYS (odd And EVEN INDEXED) WHICH REPRESENT HALF THE SAMPLING DENSITY. 
IF THE COMPARISONS ARE SIGNIFICANT* THEN ALIASING EXISTS FOR THE 
HALF-DENSITY SAMPLE. AND PERHAPS ALSO FOR THE COMPLETE SAMPLE. 

THE VALUES OF -ODDEVENS- AND -OCTAVES- SHOULD BE AS CLOSE TO 
ZERO AS NUMERICAL TRUNCATION AND/OR NOISE ALLOWS. IF THESE VALUES 
ARE GREATER THAN. SAY 1.00E-2* THEN ALIASING EXISTS FOR THIS SAMPLE 
RATE. 

A = GIVEN REAL DATA ARRAY* OF DIMENSION 2**M + 2 LOCATIONS 
S = SINE/COSINE ARRAY COMPUTED BY HARMON 
M = EXPONENT OF 2 

DATASUM - SUM OF INITIAL A ( I ) **2 FOR PARsEVALS EQUALITY 

AMPSUM - SUM OF FINAL A(I)**2 FOR PARSEVaLS EQUALITY 

ODDEVENS - ERROR BETWEEN SPECTRA OF ODD- VS. EVEN-NUMBERED POINTS 

ODDBOTHS - ERROR BETWEEN ODD-INDEXED VS ALL POINTS 

BOTHSUM - ERROR BETWEEN SPECTRA OF EVEN-INDEXED VS ALL POINTS 

OCTAVES - SUM OF OCTAVES 

DIMENSION A( 1 ) *S( 1 ) 

SET UP INDEX CONSTANTS 


N 1024 

= 2**M 



N512 = 

N1024/2 


N256 = 

N512/2 



N 102 5 

= N 1024 

+ 

1 

N 1026 

= N 1 024 

+ 

2 

N 102 8 

= N 1 026 

+ 

2 

N 1030 

= N 1028 

+ 

2 

N513 = 

N 5 1 2 + 

1 


N 5 1 4 = 

N 5 1 2 + 

2 


N516 = 

N514 + 

2 


N 5 1 8 = 

N 5 16 + 

2 


N770 = 

N514 + 

N256 

N 1544 

= N 1 030 

+ 

N 5 14 


COMPUTE THE DATA SUM FOR PARSEVAALS EQUALITY 

DATASUM = 0.0 
DO 50 I = 1.N1024 
50 DATASUM = DATASUM + A(I)*A(I) 

CALL HARM0N(A*S*M-1*-1*IFERR) 

A513 = A ( N5 1 3 ) 

A 5 1 4 = A ( N5 14 ) 

A1025 = • 5* ( A ( 1 ) —A ( 2 ) ) 

A 1026 = 0.0 


SORT COMPLEX SPECTRUM INTO ODD AND EVEN SPECTRA 
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ODDBOTH = ABS(ODD-BOTH) 

ODDBOTHS = ODDBOTHS + ODDBOTH 
EVENBOTH = ABS ( EVEN -BOTH ) 

BOTHSUM = B0TH5UM + EVENBOTH 
AMPSUM = AMPSUM + BOTH + OCTAVE 
Q = C*DC “ S*DS 
S = S*DC + C*DS 
300 C = Q 

AMPSUM = FLOAT { N512 )*< AMPSUM+A51 3*A5 1 3 + A A 1 A*A5 1 A+A1025* A 102 5*2 ♦ ) 

SORT THE SECOND HALF OF SPECTRUM 

N 1 542 = N 1544-2 
DO 600 K2 = N516»N770*2 
K 1 = K2-1 
KK2 = N 1 542 - K2 
KK1 = KK2 - 1 
AR = A ( K 1 ) 

A I = A f K2 ) 

A ( K1 ) = A ( KK1 ) 

A ( K2 ) = A t KK2 ) 

A(KK1) = AR 
600 A ( KK2 ) = A I 
C FILL IN CERTAIN LOCATIONS 

A ( N51 3 ) = A513 
A ( N5 1 A ) = ASIA 
A ( Nl 025 ) = A1025 
A ( N1026 ) = 0,0 
RETURN 
END 




Appendix D 

Subroutine DOUBLE 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


SUBROUTINE DOUBLE ( A*S*M*Nl*N2*N3) 

ROUTINE TO DO FOURIER ANALYSIS OF REAL DATA OF SIZE 2**M WHEN 
CORE CAN ONLY HANDLE ARRAYS OF SIZE 2**<M-1> 

ARGUMENT LIST — 

1. A = DATA BUFFER WORKAREA - DIMENSION ?**(M-1)+ 2 

THE EXTRA 2 LOCATIONS ARE FOR THE COMPLEX POINT AT THE 
MIDPOINT OF THE SPECTRUM 

2. S = SIN/COS TABLE - SIZE ASSUMED TO BE 2** ( M— A ) 

3. M = EXPONENT OF REAL DATA ARRAY 

4. Nl = SCRATCH TAPE WHERE INPUT DATA IS STORED IN FOUR RECORDS 

RECORD I CONTAINS A ( 1 ) ... A ( N ) 

RECORD 2 CONTAINS A(N+1) ... A(2N) 

RECORD 3 CONTAINS A(2N+1) ... A(^N) 

RECORD 4 CONTAINS A(3N+1I ... A( 4 N) 

WHERE N = 2 ** ( M— 2 ) 

5. *6. - N2 AND N3 ARE SCRATCH TAPES 

DIMENSION AU) *SI1 ) »NTAPE(3) 

NTAPEI2) = N2 
NTAPE ( 3 ) = N3 
M4096 = 2**IM-1) 

M2048 = M4096/2 
M1024 = M2048/2 
M3072 = M2048 + M1024 
M3073 = M3072 + 1 
M2049 = M2048 + 1 
M2050 = M2048 + 2 
M4098 = M4096 + 2 
MINUS = -l 
REWIND Nl 
REWIND N2 
REWIND N3 

SORT DATA INTO ODD AND EVEN ARRAYS 


DO 100 I = 1*4 

READ (Nl) ( A ( L ) * L = 1 * M2048 ) 

DO 110 J = 1.M1024 
J1 = J + M2048 
J2 = J + M3072 
A ( J 1 ) = A(2*J-1> 

110 A(J2) = A ( 2*J ) 

WRITE (N2) (AIL) *L=M2049*M3072) 

100 WRITE (N3) <A(L) »L=M3073*M4096) 

REWIND Nl 
REWIND N2 
REWIND N3 

DO FOURIER ANALYSIS OF ODD AND EVEN ARRAYS SEPARATELY 

DO 200 II = 2*3 
NT = NTAPE ( I I ) 

11 = 1 

12 = Ml 024 

DO 201 I = 1*4 

READ (NT) (A(L) * L= 1 1 * I 2 ) 

11=11+ M 1024 
201 12 = 12 + M 1024 

REWIND NT 

CALL FOURIER ( A * S *M-1 *M I NUS ) 

MINUS = -2 

WRITE (NT) ( A ( L ) *L=1*M2050) 

WRITE (NT) ( A ( L ) * L=M2049 *M4098 ) 
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non non 


200 REWIND NT 


MERGE THE EVEN AND ODD SPECTRA 


C = 1.0 

S = 0. 

DC = COS( 3*1415926536/M4096) 

DS = SIN(3.1415926536/M4096) 

DO 300 II = 1*2 

READ (N2) t At L) *L=1 *M2048 ) *ARN,AIN 

READ ( N 3 ) t A ( L ) , L=M2049 *M4096 ) * AARN * AA I N 

DO 301 K = 1 *M 1024 

K2 = 2*K 

K1 = K2 - 1 

AR = A(Kl) 

A I = A ( K2 ) 

KK1 = Kl + M2048 
KK2 = K2 + M2048 
AAR = A(KK1) 

AAI = A t KK2 ) 

DR = AAR*C + AA I *S 
D I = AAI*C - AAR*S 
AtKl) = „5*(AR+DR) 

A(K2) = .5*(AI+DI) 

AtKKl) = . 5* ( AR-DR ) 

A ( KK2 ) =• 5* t -A I+DI > 

Q = C*DC - S*DS 
S = S*DC + C*DS 

301 C = Q 

GO TO (321*322) *11 

321 A ( M2049 ) = .5 * tARN - C* ( AARN+AA I N ) ) 

A ( M2030 ) = .5 * ( -AIN + C* ( AA I N-AARN ) ) 

GO TO 323 

322 A(M20A9) = ARN + AA I N 
A ( M2050 ) = AIN - AARN 

323 CONTINUE 

WRITE (Nl) (A(L) ,L=1*M2048) 

Ml = M2 048 + 3 

MM = M4096 + M2 048 + 2 

M2 = MM/2 - 2 

DO 302 K = Ml *M2 * 2 

K 1 = MM-K 

AR = AtKl ) 

A I = A ( Kl + 1 ) 

AtKl ) = A ( K ) 

A ( K 1 +1 ) = AtK+1) 

A ( K ) = AR 

302 AtK+1) = A I 

300 WRITE ( N 1 ) (A(L) *L=M2049*M4096) 

REWIND Nl 
REWIND N2 
REWIND N3 

SORT SPECTRAL ELEMENTS INTO PROPER ORDER 


READ (Nl) ( A ( L ) * L = 
READ (Nl) ( A ( L ) * L = 
WRITE(N2> (A(L) * L = 
READ (Nl) ( A ( L ) *L= 
WRITE(N3) ( A ( L ) * L = 
READ (Nl ) ( A ( L ) > L = 

WRITE(N3) (A(L)*L= 
REWIND Nl 
REWIND N2 
REWIND N3 
READ (Nl) ( A ( L ) * L 
READ ( N 3 ) ( A ( L ) *L 

WRITE (Nl) ( A ( L ) * L 
READ ( N3 ) ( A ( L ) * L 

WRITE(NI) (A(L)*L 
READ (N2) ( A ( L ) * L 
WRITE(NI) ( A ( L ) * L 


1 9 M2048 ) 
1*M2048) 
1 *M2048 ) 
1 *M2048) 
1 * M2048 ) 
1 » M2048 ) 
1 » M2048 ) 


1*M2048) 
1 * M2048 ) 
1 *M2048 ) 
1 * M2 0 48 ) 
1*M2048 ) 
1*M2048) 
1 *M2048 ) 
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REWIND N1 
REWIND N2 
REWIND N3 
RETURN 

COMPLEX AMPLITUDES ARE ON SCRATCH TAPE Nl 
END 
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Appendix E 


Subroutine HANNING 


SUBROUTINE HANN ING ( A * N ) 
DIMENSION A(l) 

NN = N - 1 
XI = A(l) 

A ( 1 ) s XI + AI2) 

DO 100 I = 2 *NN 
X2 = A ( I ) 

All) = .5* (XI + A ( I + 1 ) ) + X2 

100 XI = X2 

A ( N ) = A ( N ) + XI 

RETURN 

END 
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Appendix F 

Subroutine HAMMING 


SUBROUTINE HAMMING! A*N) 

DIMENSION A ( 1) 

NN * N — 1 
XI - A C 1 ) 

A ( 1 ) = 1.08*X1 + . 92*A ( 2 ) 

DO 100 I =: 2 9 NN 
X2 = A( I) 

All) = .46 * ( XI + A ( I + 1 ) ) + 1.08 * X 2 

100 XI = X2 

A ( N ) = .92 * XI + 1 • 08*A ( N ) 

return 

END 
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Appendix G 


Subroutine COVAR 


SUBROUTINE COVAR (X, Y, M, S) 

C 

C X = REAL ARRAY OF DIMENSION 2**M 
C Y =REAL ARRAY OF DIMENSION 2**M 
C M =EXPONENT, 2**M =NUMBER OF POINTS 
C S =SINE/COSINE ARRAY USED BY HARMON 
C 

C CROSS-VOVARIANCE IS RETURNED IN Y 
C CROSS POWER SPECTRUM IS RETURNED IN X 
C 

DIMENSION X(l), Y(l), S(l) 

N =2**M 
N 1 =N + 1 
N 2 =2N 

DO 50 I = N1, N 2 
X fl) =0 
50 Y (0 =0 
M =M + 1 

CALL POWER (X, Y, M, S) 

DO 100 I =1, N2 
100 Y (I) = X (I) 

CALL HARMON (Y, S, M-l, 2, IFERR) 

RETURN 

END 
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Appendix H 

Subroutine POWER 


SUBROUTINE POWER ( X * Y * M * S ) 
DIMENSION X(1),Y(1)*S(1) 

N = 2**M 

CALL FOUR IER (X *S*M ) 

CALL FOURIER ( Y »S,M) 

DO 100 I = 1 * N * 2 
11 = 1+1 

TEMP = X(I)*Y(I) + X ( 1 1 ) * Y ( I 1 ) 
XIII) = X { I ) *Y (II) - X ( I 1 ) * Y ( I ) 
100 X( I ) = TEMP 
RETURN 
END 


NASA-Langley, 1968 8 
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